Shocks near Jamming 
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Non-linear sound is an extreme phenomenon typically observed in solids after violent explosions. 
But granular media are different. Right when they jam, these fragile and disordered solids exhibit 
a vanishing rigidity and sound speed, so that even tiny mechanical perturbations form supersonic 
shocks. Here, we perform simulations in which two-dimensional jammed granular packings are 
dynamically compressed, and demonstrate that the elementary excitations are strongly non-linear 
shocks, rather than ordinary phonons. We capture the full dependence of the shock speed on pressure 
and impact intensity by a surprisingly simple analytical model. 



Granular materials exhibit a wide range of complex 
collective behaviors, making them an important testing 
ground for the physics of amorphous materials (H]-fl6l]. 
The confining pressure P is perhaps the most impor- 
tant parameter controlling their properties. Strongly 
compressed granular media are, in many aspects, simple 
solids in which perturbations travel as ordinary phonons. 
However, when the confining pressure is lowered to zero, 
or the amplitude of the disturbance is much higher than 
the initial compression, the mechanical response of gran- 
ular media becomes increasingly anomalous. 

Several insights have been obtained by studying a sim- 
ple model of granular media comprised of soft frictionlcss 
spheres just above the jamming point [l[-[l6j]. The jam- 
ming point corresponds to the critical density at which 
the grains barely touch and P vanishes [l|. The first in- 
sight is that the vibrational modes of jammed packings 
resemble ordinary phonons only below a characteristic 
frequency scale uj* that vanishes as P goes to zero 
Above u>*, the modes are extended but strongly scattered 
by disorder [l3j]-[l5j]. Second, as a direct consequence of 
the nonlinear dependence of the local contact force on the 
grain deformations, the sound speed vanishes as P goes 
to zero linear sound cannot propagate when the 

particles barely touch. Third, the range of validity of 
linear response vanishes when P goes to zero. This is 
intuitive since the material is about to fall apart when 
the pressure vanishes 

As the pressure (or density) is lowered towards the 
jamming point, there are thus three anomalies that forbid 
the propagation of ordinary phonons: disorder disrupt 
phononic transport for all frequency scales, the sound 
speed vanishes and linear response is no longer valid. The 
vanishing of the sound speed and absence of a linear range 
clearly suggest that the excitations near jamming will be 
strongly nonlinear. Nevertheless, most numerical and an- 
alytical studies of energy transport have been carried out 
in solids just above the jamming point, within a vanish- 
ingly small window of linear response. By explicit design, 



these studies cannot probe non-linear energy transport 
because the dynamics of the system is solved through a 
normal mode expansion [12h15| . Therefore, with the ex- 
ception of theoretical and experimental studies on soli- 
tons in one dimensional granular chains, started with 
the seminal work of Nesterenko [l8l - l22j . non- linear en- 
ergy transport in granular packings remains largely un- 
explored. 

Numerical model. To probe how elastic energy is trans- 
ported close to the jamming point, we performed molecu- 
lar dynamics simulations of a piston-compression exper- 
iment carried out in two dimensional polydisperse amor- 
phous packings of soft frictionlcss spheres, whose radii, 
Ri, are uniformly distributed between 0.8 and 1.2 times 
their average R. Particles i and j at positions Xi and Xj 
interact via a non- linear repulsive contact potential [l2j]: 
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parameter 

effective Young's modulus of the two particles, E l} _ 
Ref. [HI for more details. The case a = 5/2 corresponds 
to Hertz's law. Lengths are measured in units of average 
particle diameter. The unit of mass is set by fixing the 
grain density to unity. The effective particle Young mod- 
ulus E* is set to one, which becomes the pressure unit. 
These choices ensure that the speed of sound inside the 
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We prepare Hertzian packings at a fixed pressure P, 
or equivalently, an average particle overlap So ~ P 2 / 3 . 
They are then continuously compressed by a piston which 
moves with a constant velocity up in the x direction 
throughout the simulation, see Fig. [T] The subsequent 
motion of the particles is obtained by numerical integra- 
tion of Newton's equations of motion subject to periodic 
boundary conditions in the y direction and a fixed bound- 
ary on the right edge of the system. We use two dimen- 
sional packings in the range of 10 3 to 10 4 particles with 
various width to length ratios. 
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FIG. 1. Snapshots of the piston-compression simulation. A 
massive piston moves to the right at a constant velocity up, 
resulting in the formation of a compression front traveling at 
a speed vs- Color indicates the local pressure at each grain. 
The average particle overlap is Sq to the right of the front and 
Ss > So to the left. 



Phenomenology. The piston compression leads to the 
formation of a front that separates two states. Ahead 
of the front, we find a region where the particles are at 
rest having the initial overlap <5o, whereas behind it there 
is a compressed region with particles moving on average 
with the piston speed up and an overlap Ss > So- Figure 
[2£i shows typical profiles for the longitudinal particle ve- 
locity u (in the x direction) as a function of x, obtained 
upon averaging velocity fluctuations in the y direction. 

Two qualitative features of the shocks stand out for all 
the amorphous packings probed in this study: the fronts 
are smooth and stable. The smoothness can be con- 
trasted with the typical shock profile that arise in ordered 
lattices of grains. Figure [2Jd, obtained for a triangular 
lattice of grains with zero initial overlap, shows large co- 
herent pressure oscillations caused by the in-phase mo- 
tion of the crystalline planes. These peaks are washed 
out by disorder in the amorphous packings. 

Second, we have systematically tested the stability 
of the front against sinusoidal perturbations (in the y- 
direction) of varying amplitudes and wavelengths in dis- 
ordered packings under various pressures. This was done 
through directs simulations [17| as well as by perform- 
ing a Dyakov's stability analysis [13, [H, [24| • A typical 
result from our simulations, illustrated in Fig. shows 
how the front remains stable due to a classic stress fo- 
cusing process, where particles "left-behind" experience 
a large compression, pushing them to catch up with the 
rest of the front. In the light of these observations, the 
shocks can be treated as one dimensional front propaga- 
tion phenomena. 



Front speed. Once transients have died out, the front 
propagates with constant speed vs in the amorphous 
packings. Upon using conservation of mass across the 
shock front, we derive a one dimensional relation be- 
tween the characteristic velocities up and vs, through 
the average radius of the particles, R, and the average 



v s = Up 



2R-S 
Ss - S 



(2) 



Since the particle compression Ss is typically much less 
than its diameter 2R, Eq. ([2]) implies that vs ^ up. This 
is consistent with our numerical findings summarized in 
Fig. [3^ where the dependence of vs on up is explored 
systematically for different compressions. 

Inspection of Fig. [3^ reveals two distinct regimes. Fol- 
low up, the front speed vs is nearly independent of up - 
in this (quasi)linear regime, vs is simply controlled by the 
initial pressure P. The strongly non-linear shock waves 
regime is reached for high compression speed up, where 
vs depends on up, but not on P. 

The data for vs can be collapsed onto a single mas- 
ter curve, as shown in Fig. [3Jd. We achieve this upon 
rescaling the vs axis by fs(0), the numerically deter- 
mined value that the front speed attains in the limit 
of vanishing up (sec Fig. [3^). The up axis is rcscalcd 
by a pressure-dependent velocity scale up, obtained by 
matching the low and high up asymptotes in Fig.[3K (see 
arrow): up marks the crossover between linear acoustic 
waves and shocks. 

Scaling analysis. The pressure dependence of vs(0) can 
be rationalized using scaling arguments. We expect that 
vs(0) reduces to c, the speed of linear longitudinal sound 
waves. To determine the scaling of c with pressure, note 
that c ~ vfi, where the bulk modulus B = ^ and P = 
The change in volume dV scales linearly with Sq, 




FIG. 2. (a) Profiles of a shock wave at different times ob- 
tained by averaging the particle velocity in the y direction 
(symbols). The red lines shows the fits of the fronts to an em- 
pirical fit formula, (b) Oscillatory velocity profile of a shock- 
like wave propagating through an hexagonal array. The inset 
shows the in-phase motion of the crystalline planes that gives 
rise to the pressure oscillations, (c) Representative snapshots 
of the focusing and flattening of an initially curved front gen- 
erated by a sinusoidal piston (time progresses from left to 
right). 



the average overlap between particles, while the energy 
scales as E ~ Sq , see Eq. Q] Upon setting a = 5/2, we 
obtain the pressure dependence of the longitudinal speed 



of sound c 



,1/4 



P 1 / 6 valid for Hertzian interactions 



[121 ] . Figure [3j3 shows that the numerical data for vs(0), 
represented by red symbols, is consistent with the S^ 4 
scaling, which is shown as a continuous red line. 

We now turn to the regime of high piston speeds, 
tip > Up, when the front speed v$ becomes nearly in- 
dependent of P. Since up, R and (5 are all known, we 
need one additional relation which combined with Eq. ([2]) 
will make a definite prediction for the shock speed. We 
note that for strong shocks, the propagating front gen- 
erates a characteristic compression S ^> So and a cor- 
responding increase in the kinetic energy. By assuming 
that the kinetic and potential energies are of the same 
order, we obtain up ~ J 5 / 2 . We have tested numeri- 
cally that this non-trivial proportionality relation exists 
for strong deformations, see Fig. [5J1. Upon combining 
the balance between kinetic and potential energy with 
Eq. @, one readily obtains the power law v$ 



i 1/5 

ip , 



plotted as a dashed line in Fig. \Sjp. This scaling rela- 
tion is clearly consistent with our numerical data for the 
speed of strongly non-linear shock waves. 

We deduce the dependence on compression of the 
crossover speed u* P by smoothly matching the two asymp- 
totic relations for the front speed vs ~ w)/ 5 and i>s(0) 



,1/4 



This leads to the power law relation up 



,5/4 



(continuous blue line in Fig. [3J;) that is consistent with 
our numerical values (blue symbols). Note that the 
data collapse in Fig. [5b depends only on the scaling 
up ~ (j[j/ 4 and it is not sensitive on the precise defini- 
tion of the crossover speed. Upon using the conversion 
relation up ~ <5 5 / 4 , the intuitive expectation that the 
crossover takes place when S ~ <5o is confirmed. 

We conclude that by controlling So or P, which param- 
eterize the distance to the jamming point (at P = and 
So =0), we can tune up and the onset of the strongly 
non-linear response of the packings. Our key numerical 
findings on the shock velocity summarized in Fig. [3J can 
be grasped from scaling near the jamming point. 
Analytical model. In order to account for the depen- 
dence of vs on up and the smoothness of the shock pro- 
files, we construct the simplest possible ID model that 
quantitatively accounts for the trends observed in Fig. [3] 
and sheds light on the role of disorder. 

In the continuum limit, we obtain the following equa- 
tion governing the dynamics of the system in terms of 
the strain field 5(x,t) [25| : 



^ ®ttxx 



Stt 



4R 2 e 



[s a -% x = 0. 



(3) 



To gain some intuition for the physics behind Eq. 
note that by setting a = 2, one recovers a linear disper- 
sive wave equation, with speed proportional to y/e/m in 
the long wavelength limit. By contrast, when a > 2 a 
non-linear wave equation is obtained. Nonlinearities and 
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FIG. 3. a) Speed of the front vs versus particle velocity up 
measured in units of v g , the sound speed within the grain, for 
decreasing particle overlap So. b) Same plot as in (a) but with 
vs normalized by vs(0) and up normalized by the crossover 
particle speed up: fs(0) and up are indicated in panel (a). 



The dashed line indicates the power law vs 



Up 



l/s 



charac- 



teristic of a sonic vacuum. The black line indicates the the- 
ory developed here to describe the universal transition from 
weakly to strongly non-linear waves in systems close to jam- 
ming, c) Variation of i>s(0) and up with distance to the jam- 
ming transition parameterized by the initial average overlap 
5q. The dashed lines indicate the power laws vs(0) ~ 5$ , 
up ~ 8f/ 4 . d) Variation of the kinetic energy with potential 
energy in dimensionless units - same color code as in (a-b). 
The dashed line indicates the linear relationship observed for 
strong shocks. 



dispersive effects gives rise to finite amplitude waves: ei- 
ther solitary waves or shocks are possible depending on 
the drive [19J. 

Shock propagation is modeled by the combined strain 
S(x,t) = So + g(x), where g(x) gives the shape of the 
shock and x = x — v$ t. Upon inserting this ansatz into 
Eq. ([3]), we obtain the conservation law i<5| + W{8) = 0, 
where W{8) is given by 



W{5) 



24e 



motv 



S 

24<V 



^% 5 o- 2 -^)(S-So). 



(4) 



This conservation law can be interpreted as describing 
the total energy of an effective particle at position S 
rolling down a potential well W(5), shown as a red line 
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FIG. 4. (a) Simulations of an ordered chain with small vis- 
cosity display large coherent oscillations in the front profile 
(black line) [20]. If the viscosity is large enough, one obtains 
an homogeneous shock profile, shown as a green line, similar 
to the profile in Fig. [2^. (b) The presence of an effective 
viscosity will induce the oscillation of the particle (black tra- 
jectory) towards the bottom of the potential W(&), shown as a 
red line. If the viscosity is large enough, the particle can move 
directly to the minimum without performing any oscillations, 
see the green trajectory corresponding to the homogeneous 
shock profile of panel (a). 



in Fig. [4£i (here x maps to time so that \b\ is the kinetic 
term of the particle) [2a |. 

One of the key ideas of our work is that disorder can 
act as an effective viscosity for the shock: the energy 
imparted unidirectionally by the piston is redistributed 
among other degrees of freedom, reducing the energy 
propagating with the shock front. In our mapping, this 
implies that the effective particle, initially located at the 
maximum of the potential W = 0, moves to the minimum 
of the potential well (see Fig. 0k) . Thus, upon setting 
dsW(S) = 0, we can obtain a relation between propaga- 
tion velocity and induced compression in the front 



1 (Ss/So)*- 1 - 1 



Q 



1 (Ss/S ) - 1 



(5) 



that is independent of viscosity, even if an infinitesimal 
amount of dissipation is necessary to obtain a steady 
state solution of Eq. (|3]). 

Together, Eqs. (J5|) and §5§ can be seen as a parametric 
relation between front and particles velocities, where the 
overlap 5s produced by the passage of the front is the 
parameter. Such a parametric plot of vg versus up is 



drawn as a continuous curve on the numerical data in 
Fig. This comparison shows that Eqs. © and ([5]) are 
in excellent agreement (without any fitting parameter) 
with the results of our numerical experiments on shock 
propagation. 

Discussion. The shock formation explored in the 
present study is a generic phenomenon independent of 
the dimensionality of the sample that relies purely on 
the presence of a nonlinear law between grains (for any 
a > 2) and not on the presence of friction. Experimen- 
tally, this can be tested by impacting a box of (frictional) 
glass beads with a heavy mass, for a range of impact 
speeds and pressures - preliminary experimental results 
for the front speed compare favorably to our theoretical 
predictions in Fig. l3l [271]. 

We note, however, that in frictional granular media, a 
second type of densification front can be observed, which 
is often referred to as plowing [i| ■ Whereas our shock 
waves always propagate with speeds above the linear 
sound speed, and continue to propagate even after the 
driving stops, plowing fronts are generally much slower 
(in of order 1 m/s), and stop almost immediately 
when the driving stops. We believe that the underlying 
difference is that our shocks are dynamical phenomena, 
set by a balance of potential and kinetic energies, whereas 
plowing is in essence a quasistatic phenomenon, domi- 
nated by dissipation. In the dynamic case, the change 
in packing fraction induced by the shock is associated 
with grain deformations, whereas in the quasistatic case, 
densification is dominated by grain rearrangements and 
compaction. 

Outlook. The shocks that arise in grains near 
jamming are just one representative of a broader class 
of strongly nonlinear excitations that emerge near the 
marginal state of suspensions, emulsions, wet foams and 
weakly connected fiber networks @, [28|, [2{|. Close to 
losing their rigidity, all these materials exhibit a vanish- 
ing range of linear response, so that almost any amount of 
finite driving will elicit an extreme mechanical response 
in the form of rearrangements, yielding and flow [l6j], 
[30j ] . j3l| . [32|]. It remains an open question whether all 
these phenomena can be successfully described in terms 
of simple scaling near jamming. 
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